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Abstract 

We present a novel treatment of finite temperature properties of the one- 
dimensional Hubbard model. Our approach is based on a Trotter-Suzuki mapping 
utilizing Shastry's classical model and a subsequent investigation of the quantum 
transfer matrix. We derive non-linear integral equations for three auxiliary functions 
which have a clear physical interpretation of elementary excitations of spin type and 
charge excitations in lower and upper Hubbard bands. This allows for a transparent 
analytical study of certain limiting cases as well as for precise numerical investiga- 
tions. We present data for the specific heat, magnetic and charge susceptibilities for 
various particle densities and coupling strengths U. The structure exposed by these 
curves is discussed in terms of the elementary charge and spin excitations. Special 
emphasis is placed on the study of the low-temperature behavior within our ab initio 
approach confirming the scaling predictions by Tomonaga-Luttinger liquid theory. 
In addition we make contact with the "dressed energy" formalism established for 
the analysis of ground state properties. 



*e-mail: gj@thp.uni-koeln.de 

' e-mail : kluemper @ thp . uni-koeln . de 

■^e-mail: suz@hepl.c.u-tokyo. ac.jp, Permanent address: Institute of Physics, University of Tokyo at 
Komaba 



1 



1 Introduction 



The Hubbard model represents the most fundamental model for highly correlated electron 
systems. It has therefore attracted much attention since its formulation. Triggered by the 
discovery of high T c superconductivity, correlation effects of the model have been of strong 
recent interest. 

For the problem in one spatial dimension an exact solution is available via the Bethe 
ansatz It clarified the Mott insulator nature of the itinerant electron system at half- 
filling in ID. The extensive list of investigations of zero temperature properties ranges 
from studies of the elementary excitations [@, || £|, [|, over magnetic properties || to 
studies of correlation functions in the strong coupling limit J?|. Recently, the large group 
of symmetries of the Hubbard model was analyzed in || |], [K], [LI]]. It was, however, almost 
20 years after the discovery of the exact solution that various properties like asymptotics 
of correlations functions at T = have been investigated by using results of conformal field 
theory. The Tomonaga-Luttinger liquid properties of the one- dimensional model have been 
shown not only at a qualitative but also at satisfactory quantitative level |12|, [13], [L4], [15|] . 

In this report, we address the problem of the Hubbard model at finite temperatures. 
In fact, soon after the seminal solution a thermodynamic formulation was set up 
through the string hypothesis approach [16]. It results into infinitely many (oo x oo) 
coupled nonlinear integral equations for infinitely many unknown functions. Obviously, 
quantitative studies of such equations need much effort and were performed only relatively 
recently (also about 20 years after the thermodynamical formulation!) [[17] Still, the 
explicit calculation allowed for only 2 (!) bound charge rapidities and 15-30 spin rapidities, 
while the original equations contain oo x oo rapidities. 

Here we attack the problem via a completely different approach developed recently 
which avoids the computational complications and also renders correlations lengths at 
finite temperatures accessible. We make use of a general equivalence theorem between 
(^-dimensional quantum systems at finite temperatures and d + 1-dimensional classical 
systems [ID], thus we employ a convenient mapping to a two-dimensional classical model. 
The evaluation of the free energy then reduces to finding the largest eigenvalue of the so- 
called quantum transfer matrix (QTM) |g [2g, |2T], ^, |23], [?§ g§ g§ £7], |2§ . The crucial 
observation is the existence of a commuting family of QTMs labeled by one complex 



parameter (spectral parameter) |25|, |26], |3(J, plf| . This makes the meaning of integrability 



manifest, and allows for the investigation of thermodynamics through the study of the 
analytical properties of suitably defined auxiliary functions on the complex space of the 
spectral parameter. One of the most practical advantages in this novel formulation is the 
fact that one has to deal with only a finite number of auxiliary functions and nonlinear 
integral equations among them. Therefore, we can expect results with higher numerical 
precision. Furthermore, the involved auxiliary functions have clear physical interpretations 
in terms of elementary excitations at T = 0. 

Such a strategy has been successfully applied to several interesting models including 
the spin 1/2 XYZ model and derivates EE EM EM E7J, the integrable t - J model P, 



the super symmetric U model [31], and (with limited success) to the Hubbard model 
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32| , |23| , The main restriction of the previous work on the Hubbard chain |23|, Q 
is the limitation to the case of half-filling. For any finite doping the resultant equations 
appeared to be numerically ill-posed. Here we want to overcome the technical difficulties 
and derive a set of equations that allow for convenient (numerical) studies and clear 
analytical insight for all particle densities. Before doing so, we want to remind of the 
remarkable differences in comparison to other solvable models. The R— matrix for the 
classical analogue ( "Shastry's model" ) pi), |35], does not possess the difference property 
of rapidities. Therefore, the intertwiner depends on two spectral parameters, not only on 
the difference. Such violations are only known for Shastry's model and the chiral Potts 
model. In view of analyticity, the Hubbard model is also quite unique. One can easily 
recognize this by comparing the Bethe ansatz equations (BAE) at T = for the Hubbard 



model PJ, with, for instance, those for the integrable t — J model p?|, |38], |39|. In both cases 



there are two kinds of BAE roots corresponding to charge and spin degrees of freedom. 
However for the integrable t — J model both types of roots vary from — oo to oo, while the 
charge-rapidities for the Hubbard model only vary from — it to tt with a corresponding 
periodicity. This different character of the BAE roots brings about branch cuts in a 
complex parameter plane as we will see below. The roots show an exotic behavior: they 
flow from one Riemann sheet to the other with changing temperature. 

Despite these difficulties, we will show our strategy is successfully applicable to the 
Hubbard model (and finally overcomes the technical problems still left in the earlier 
approach With a careful choice of auxiliary functions, we can completely encode 

the information about zeros in both sheets. They are shown to have close relation to the 
physical excitations of holons and spinons in the T — > limit. This limit will be studied 
in quite some detail as it allows for a first principles derivation of Tomonaga-Luttinger 
liquid properties at low but finite temperatures. Several quantities of physical interests 
are evaluated with high numerical precisions for wide ranges of temperatures and fillings. 

This paper is organized as follows. In Section 2 we present the mapping of the Hubbard 
chain at finite temperature to a two-dimensional classical system with integrable QTM. 
In Section 3 the eigenvalue equations for the QTM are derived which are cast into a 
difference type form in Section 4. Section 5 is devoted to the derivation of non-linear 
integral equations for Hubbard interactions U > 0. In Section 6 the integral equations are 
investigated numerically and the results are discussed. Section 7 deals with the analytical 
study of various limiting cases, notably the low-temperature asymptotics. In Section 8 
we present our summary and outlook. The derivation of the integral form of the QTM 
eigenvalue is deferred to Appendix A. 



2 Shastry's Model as a classical analogue of the ID 
Hubbard Model 

In the novel formalism, it is essential to deal with the two dimensional classical counter- 
part. Fortunately, Shastry has already found two classical versions for the Hubbard model 
f40| , |34|. For the latter model, a proof of the Yang-Baxter integrability has been given in 
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a recent analysis pTJ by the use of the tetrahedron algebra. As the Yang-Baxter relation 



makes the finite temperature analysis much easier, we adopt the latter version here. 

We sketch the essential properties of the model. The Hubbard model describes a 
lattice fermion system with electron hopping term and on-site Coulomb repulsion with 
Hamiltonian 

^Hubbard = \ -( C l+l,v C i,« + 4,^+1,^) + U ("i - ~ ~ \) 

1=1 {a=± 

"^external • ( 1 ) 

The external field term H ex ternai = — SJM n «,++ n j ,-)+H/2(n it+ — rii _)] will be omitted for 
the time being. According to ||34j| , it is easier to find a classical analogue after performing 
the Jordan- Wigner transformations for electrons in ID. The resultant spin Hamiltonian 
is 



n,n+l 5 
n=l 

7~in,n+l — ( a n a n+l + a n+l a n ) + { T n T n+l + T n+l T n ) + T a 'n T ni (2) 

where L denotes the chain length of the system. Note that we are now imposing peri- 
odic boundary conditions for the spin system (o"i(r 1 ) = a ■ This does not give 
the periodic boundary conditions for the underlying electron system. The differences in 
boundary conditions, however, will not affect thermodynamic quantities like the specific 
heat. 

For the counterpart in two dimensions, we consider two double-layer square lattices, 
say a o and a r lattice. Each edge possesses a local variable ± and each vertex satisfies 
the ice rule. The vertex weights consist of contributions from both intra and inter lattice 
interactions. The intra part is given by the product of vertex weights of the free-fermion 
six vertex model: £i i2 [u] = ii 2 (u) ® l\ 2 (u) where 

nu , \ a ( U ) + K U ) a ( U ) - K U ) z z I + - - +\ /on 

*i, 2 («) = 2 2 12 2 + Gl ^ ® 

and a(u) = cos(m), b(u) = sin(u), c(u) = 1. Taking account of inter- layer interactions, the 
following local vertex weight operator (denoted by S) is found J34J 



Si,2(v,u) = cos(u + v) cosh(h(v,U) — h(u,U)) £ 1)2 (v — u) 

+ cos(i; — u) sinh(/i(t>, U) — h(u, U)) i\^{u + v) a^r^ 

where sinh 2h(u, U) := Ua(u)b(u)/2. The Yang-Baxter relation for triple S matrices is 
proved in [^T|. The commutativity of the row-to-row transfer matrix, 

T(u):=f[S^(u } 0) (4) 



4 



is the direct consequence. 

The S matrix and 7i are related by the expansion in small spectral parameters, 

S lj2 {u, 0) = S 1>2 {0, -u) ~ P(l + un 1>2 + 0(u 2 )), (5) 

where P denotes the permutation operator, P(x®y) = y®x. Once the Yang-Baxter rela- 
tion is established, we can apply our machinery for thermodynamics. Here we summarize 
the necessary formulas, see |3(J for details. We define Ri t2 (u, v) = Si t2 (v,u)\u-,-u, and in- 
troduce R(u, v) and R(u, v) by clockwise and anticlockwise 90° rotations of Rx j2 (u, v). We 
further introduce an auxiliary transfer matrix T(u) made of Boltzmann weights -R(0, — u). 
The partition function is given by 

Z = lim Tre"^ = lim lim Tt [T(u)T{u)] N/2 \ u=m (6) 

where H'l differs from TCl by the sublattice gauge transformation, c„ )CT — > ( — l) n c„, jCr which 
does not affect thermodynamic behaviors. We regard the resulting system as a fictitious 
two-dimensional model on a L x iV square lattice, where N is the extension in the fictitious 
(imaginary time) direction, sometimes referred to as the Trotter number. Now by looking 
at the system in a 90° rotated frame, it is natural to define the "quantum transfer matrix" 
(QTM) by 

N/2 

T qTyi (u,v)=<S)R(-u,v)®R(v,u). (7) 



The interchangeability of the two limits (L, N — > oo) leads to the following expression, 

1 L 



lim lim Tr 

N— >oo L— >oo 



There is a gap between the largest and the second largest eigenvalues of 7q TM (m,0) for 
finite f3. Therefore, we have a formula for the free energy per site, 

/ = -k B T lim In A max (u = ^ ) . (9) 

where A max (w,0) denotes the largest eigenvalue of 7q TM («,0). Now the evaluation of 
the free energy reduces to that of the single eigenvalue A max . Of course, a sophisticated 
treatment is necessary in taking the Trotter limit —>■ oo as u now explicitly depends on 
it. The following sections are devoted to this analysis. 

A general comment is in order. It seems somewhat redundant to define 7q TM (u, t>) as 
we only need the value at v = 0. This formulation, however, manifests the integrability 
structure and the existence of infinitely many conserved quantities. This is best seen in 
the commutativity of transfer matrices 

Pqtm(«, v), TqtmO, v')] = 0, (10) 

with fixed u. One can prove this by showing that two QTMs are intertwined by the same R 
operator as for the row-to-row case. The outline of the proof is graphically demonstrated 
in Fig.|l]. The existence of the parameter labeling the family of commuting matrices makes 
the subsequent analysis much more transparent. 
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Figure 1: a) Graphical depiction of the fundamental Yang-Baxter equation for R and the 
associated one for R and R obtained through rotation, b) "Railroad proof" for the com- 
mutation of two QTM's with any spectral parameters v and v'. Due to a) the intertwiner 
for R vertices is identical to the intertwiner for R vertices. 



3 Diagonalization of the Quantum Transfer Matrix 

It has been an issue of current interest to find an explicit algorithm of the diagonalization 
of the row-to-row transfer matrix of Shastry's model or its fermion version. The standard 
tool in such studies, the quantum inverse scattering method, has been applied and turned 



out to be successful after elaborate calculations J42|. We also refer to the analytic Bethe 
ansatz study j|3 . 



Here we want to diagonalize the QTM. At first glance, the diagonalization scheme 
for this looks quite different from the row-to-row case. The QTM has a complicated 
inhomogeneous structure seemingly demanding much more effort. Fortunately, this is not 
true. The crucial observation is, as remarked in the previous section, that QTMs share 
the same intertwining operator with the row-to-row transfer matrices. In view of QISM, 
this results into identical operator algebras which allows for the diagonalization of the 
trace of the monodromy matrix. (Note that we adopt periodic boundaries in the Trotter 
direction.) Thus, the eigenvalue equation of the QTM involves the same combinations 
of "dress functions" (in the terminology of the analytic Bethe ansatz) as in the row-to- 
row case. One only has to replace the vacuum expectation values taking account of the 
quantum space and the inhomogeneity. 
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We define state vectors i — 1, • • • , 4 by 

|1> = !+,->, |2> = |+,+>, |3> = |-->, |4> = |-+>. (11) 

A convenient vacuum in the present study may be \Q) := |1,4, 1,4, ■ ■ • , 1,4). Then the 
vacuum expectation values read 

A = (Q\T ht \Q) = {^{-u,v)Bi%v,u)) N ' 2 i = 1,... ,4. (12) 

Substituting explicit elements for R, we find the relations 



Ai/A 



2 



N/2 



A 4 /A 2 
A 2 



(1 - z-{w)z+(x)){l - z+(w)z + (x)) \ 
(1 + z4w)z + (x))(l + z + (w)z + (x))J 
A 2 = A 3 

(1 + z-{w)/z-{x)){l + z + {w)/z-{x)) \ N/2 



(1 — Z-(w) / z-.(x))(l — z+(w) / Z-(x)) 



(1 + e ix ){l + e 4 ™) v ^_(w) Z-(^) ^ ' z_(w)' 

where we have introduced the parameterizations x, w for v, u, 

e 2x = tanv, e 2w = tan it, (13) 

and the functions 

z± (x) = e 2M-)±2x ; 2h ^ = _ sinh -i ( U \ _ (M) 

\4cosh2x / 

Now that we have the explicit vacuum expectation values, the eigenvalue can be written 
down directly thanks to the above argument 

A ( W ) _ J3(*+H/2) A l TT 1 + ZjZ-[x) 



TT 2a; 

a 2 ^ 2 ii 6 r-^z + (x) 



, 2^ TT 2x 1 + z 3 z -( x ) TT M^) - l/z_(x) - 2iw a + 3U/2 
1 = 1 1 - ZjZ+(x) 1 = 1 - l/z-ix) - 2iw Q + U/2 

2x 1 + ^+ (x) / 2j TT- z_ (x) — 1/ z_ (x) — 2iu> Q , — U/2 



(15) 



TT ^ TT 

^ 11 i _ „ M/,. 11 



z_ (x) / Zj ,2_ (x) — 1/ (x) — 2iw a + U/2 



■ 1 V" / / *7 1 

i 0(n-H/2) ^±Y\ p~ 2x 1 + Z +( X )/ Z i 



Note that we imposed a non-vanishing chemical potential and magnetic field at the last 
stage, as they merely lead to trivial modifications in A due to twisted boundary conditions 
for the QTM g(| . 
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The parameters {zj}, {w a } satisfy the BAE, 



e P(n,-H/2) 



;i + ^_H/^)(i + ^ + H/^)\ iV/2 _ , 1W A z j -i/ Zj -2iw a -u/2 



n 



e 



z-(w)/z s ){l-z+(w)/zj)) ^ Zj-l/zj-2iw a + U/2- 



. =1 ^ - l/zj - 2iw a - tf/2 ~ ~ 11 2i(w a -wp)+U' 

(16) 

Here some remarks are in order: 

1. Although the validity of expression ( |l5f) is a logical consequence, it would be a good 
exercise to verify this form for one-particle states. For example, we take 7^ 4 (i/)|f2) as 
a representative, and calculate its eigenvalue. A standard argument in QISM leads 
to the following "wanted terms" . 



wanted terms 



, RsA v ^) . . ^4,4 
+ 3 ^37 — r + A 4^37 — r 



(17) 



^,4(^)|n> 



A straightforward however lengthy calculation shows the coefficient in ([T7p is equal 
to (P|) with m = 1, £ = 0, z x = z_(l/21og(tauz/)). 

2. We have verified that fll5[) gives the largest eigenvalue identical to the one obtained 
by brute force diagonalizations of finite systems up to A = 6. The groundstate lies 
in the sector m = N, £ = N/2. For the repulsive case and \i = H = 0, z/s are all 
on the imaginary axis, while w^s are real. 

3. The free-fermion partition function is easily recovered for U — 0. 

4. Starting from another vacuum \Q') = |2,3, ■ ■ ■), one reaches a different expression. 
The resultant one is actually equivalent after negating U and exchanging H/2 <-> fi, 
namely, after the partial particle-hole transformation. 



4 Associated auxiliary problem of difference type 

The thermodynamics leading to the free energy is encoded in the solution to the BAE 
in the limit N —>■ 00. For finite N it is possible to solve the BAE numerically. However, for 
large N it is quite complicated to find the numerical solution even for the ground state. 
Furthermore, in the Trotter limit N — > 00 the roots {v^,w^} accumulate at infinity. This 
is similar to other models (Heisenberg model, t — J model) where the solutions of BAE 



of the QTM tend to the origin pH 061 27, 29, 50, 31]. It represents the main problem in 
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analyzing the limit N — > oo directly on the basis of the BAE. To overcome this difficulty 
one can express the solution of the BAE by a system of non-linear integral equations. 
This has been done for several models E3, H5|, p5L 0, 53, 130, |3_1]. 



The first problem to be overcome is the involved structure of BAE (0). Upon intro- 
ducing the quantities 

1 



Sj ~ 2i \ Zj z 



the equations fllq) take a difference form in the rapidities {sA, {w a }: 



with 



and 



e-^-W^) = - q2 [ S \% (19) 
q 2 {w a - 2ij) qi(w a - ij) ' 



Qi{s) = Y[{s- Sj), q 2 (s) = JJ(s - w a ), 7 = -, (21) 



(i-^MM 8 ))(i- z+ M/z(»)) \ y/2 
+ z (#W)( ihWM«))J ' 1 ' 

z(s) = is(l + Va-l/* 2 )). (23) 

The function z(s) possesses two branches. The standard ("first") branch is chosen by the 
requirement z(s) ~ 2is for large values of s, and the branch cut line [—1, 1] (corresponding 
to a cut for ^fz from — oo to 0.) We will not refer explicitly to the second branch of z(s) in 
this work. However, for <j)(s) both branches will be used. In order to distinguish between 
the first and second one we will use the notation <p + (s) and <f)~(s), respectively. 

In Fig.|2] the distribution of rapidities Sj is shown for iV — > oo. Note that there are 
infinitely many rapidities on the first (upper) sheet and finitely many on the second 
(lower) sheet. The number of rapidities on the second sheet is increasing with decreasing 
temperature thus resulting into a flow from the first to the second sheet. 

Note that the general expression (|15|) is quite complicated, but simplifies considerably 
at v — and u — > 



in 



A(v = 0) = e^(l + e ^ +H ^)(l + e^~ H ^)u N ]J Zj. (24) 

3=1 

Lastly, we want to comment on the difference type property of (0). These equations 
are "BAE compatible" with the following eigenvalue equation of an "auxiliary transfer 
matrix" 

A aux ( S ) = Ai(s) + A 2 (s) + A 3 (s) + A 4 (s), (25) 
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Figure 2: Depiction of the flow of rapidities Sj from the first (upper) Riemann sheet to 
the second (lower) one for decreasing temperature. 



Ai(s) = e^m^^l , A 2 ( S ) = e ^ " 2 [ S ~ 2l7 \ , 
Hs) = fP^- , X,(s) = l - 



<h{s)(h{s + ilY " 4>(s + i"f) qi(s + ry) 



The reason is obvious, A aux (s) is pole free under BAE (|20|) just like the original QTM. 
The construction (|26| ) at this point is purely mathematical, however it will be the starting 
point of the derivation of integral equations in the next section. 



5 Non-linear integral equations 

In this section we are concerned with the derivation of well posed integral equations equiv- 
alent to the nested BAE for the largest eigenvalue of the QTM. (We restrict ourselves to 
the case U > and point out that U < is simply obtained via a particle-hole transfor- 
mation.) We introduce a set of auxiliary functions which satisfy a set of closed functional 
equations which later on are transformed into integral form. The explicit expression of 
the functions b, b, c, c which proved useful is 

. h + h + h + k j- h + k + h + k 

b = : : — , b 



h + k + h + k h + h + h + h 
_ h + k h + k + h + h 

^3 + ^4 /1 + /2 + /3 + ^4 + ^1+^2+13+ IT 

_ _ h + h h + h + h + k 



(27) 



h +l 2 h + h + h + h + h + h + h + k'' 
where the functions lj and lj are closely related to the Xj defined in ( 26) 

fj( a ) = A i («-i 7 )-^VW(a), (2g) 
lj(s) = Xj(s + 27). 
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The main observation in connection with the functions defined in fl2"T|) is based on 
elementary facts of the theory of complex functions. In particular any analytic function 
on the complex plane is entirely determined by its singularities, i.e. poles and branch cuts, 
as well as its asymptotic behavior at infinity. Below we will show that the singularities of 
In b, In c etc. on the entire complex plane are exhausted by the singularities of ln(l + b), 
ln(l + c) etc. close to the real axi^. Furthermore all the involved functions show constant 
asymptotics for N finite. Hence there exists a suitable integral representation of In b, In c 
etc. in terms of ln(l + b), ln(l + c) etc. The latter functions will be abbreviated by 

m t , . h + h + h + k + h + h + h + U 
23 = 1 + b = 



23 = 1 + b 
€= 1 + c 

£ = 1 + C : 



k + k + k + k 
k + k + k + k + k + k + k + k 

k+k+k+k 
k+k+k+k k+k+k+k+k+k 



(29) 



k + k k + k + k + k + k + k + k + k' 
k+k+k+k k+k+k+k+k+k 



k +k k+k+k+k+k+k+k+ k 

Quite generally all the above auxiliary functions have a product representation with factors 
of the type ... + l 3 + k + k + k + ■■■ • As a matter of the BAE the poles of each lj and 
lj function in ... + k + k + k + k + ■■■ are canceled by the neighboring terms. Poles can 
only "survive" if such a string does not begin with l\ or does not end with k- There are 
extended singularities (cuts) due to the function <fi appearing in the definition of Ai and 
A4. Hence all terms k + k + ■■■ and ... + k + k possess branch cuts along [—1,1]+ 2ry and 
[—1, 1] — 2«7, respectively. Furthermore, terms like ... + k + k and k + k + ■■■ have branch 
cuts along [—1,1]. However in combinations ... + k + k + ■■■ the branch cut due to the 
function disappears, because 

k(s) + k(s) = S^ ^±^ , (30) 

and <f) + {s) + <f>~(s) is analytic everywhere. 

Inspecting the function Ai + A 2 + A 3 + A 4 more closely we find poles of order N/2 at 
s — i-f and 27 — s where 

z(s ) = z-(w), 2iso ~ Nj 13 for large N. (31) 

In addition we find zeros and branch cuts on the lines Im(s) = ±7 which we denote by 

N 

ln[Ai(s) + A 2 (s) + A 3 (s) + A 4 (s)] =a-~2 ln [( s + *7 - s )(s + s - 27)] 

+ L_(s + i r i) + L + (s - 27), 



1 The relevant singularities are distributed exactly on the real axis for vanishing external fields. For 
this case the subsequent treatment can be taken literally. For finite external fields h, [i deviations from 
the real axis occur. The following reasoning still applies mutatis mutandis. 
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where = s denotes that left and right hand sides have the same singularities on the entire 
plane, and L± are suitable functions possessing the desired singularities and being ana- 
lytic otherwise. (The existence of such functions can be proved quite easily. An explicit 
expression is given by contour integrals of the type fl37D -) From (|32|) we find the following 
singularities 



Hhis) + l 2 (s) + l 3 (s) + k(s)} = s - — ln[(a - s )(s + s - 2z 7 )] + ln[0 + (s)<T(s)] 

+ L_(s) + L + (s-2i 7 ), 

— — — — N 
ln[^(s) + l 2 (s) + l 3 (s) + k(s)} =s~— ln[(s + 2i 7 - s )(s + s )} 

+ L_(s + 2i 7 ) + L+(s), 
From this, and fl2"7| , |2D| ) and the identity 

(a - s ){s + s - 217) 



(33) 



_(s + s )(s - s + 2z 7 )_ 



N/2 



(34) 



we find the singularities 



In b(s) = s L_(s + 2ry) + L + (s) - L_(s) - L + (s - 2i 7 ), 
lnQS(s) = s - L_(s) +rest, (35) 
lnc(s) — lnC(s) = s L_(s) — L+(s) + rest, 

where "rest" indicates singularities not located on the real axis. 
Next we introduce the notation 

(gof)(s) = J g(s-t)f(t)dt, (36) 

for the convolution of two functions g and / with contour £ surrounding the real axis at 
infinitesimal distance above and below in anticlockwise manner. From Cauchy's theorem 
we find for any function / analytic above and below the real axis 

(k o f)(x ± ie) — (k o f)(x) + f(x ± ie), where k(s) = — ;-. (37) 



2iri s 



For further convenience we introduce the functions 

7 /vr 



Ktis) = k(s) -k(s + 2i-f) 

ays - 

Kl(s) = -k(s) + k(s - 2* 7 ) = l/7 [. (31 

S{S — 2vy) 



s(s + 2i 7 ) 
V 

K 2 (s) = k(s - 2t 1 ) - k(s + 2 %1 ) = -j 2 ^ 
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which will play the role of integral kernels. From (|35|J37|J38| ) we find 

[K 2 o In 03 + K[ o (lnc - InC)] = s + 217) + L+(s) - L_(s) - L + {s - 2i 7 ). (39) 

Upon comparing (^,^) we conclude 

lnb(s) = if 2 °ln03 + ^o(lnc-in(_) + const, (40) 

as both sides are complex functions with identical singularities. (The difference function 
is entire, i.e. analytic on the entire complex plane. Furthermore the difference function is 
bounded, hence it is constant.) The constant is computed from considering the asymptotic 
behavior at s — > 00 

const = -/3H. (41) 

For the derivation of the second type of integral equation we define an intermediate 
set of auxiliary functions 

h + h _ 1 h+h+h+h 

r = rrT' t = i + t = — — — — , 

_ _ ___ ___ (42) 

_ h + h 7p . , _ h + h + h + h 
r== — =, T=l+r= = — = . 

h +h h+ h 

Quite similar to the above derivation we conclude the identity 

lnr(s) - — In — -^_2 — + pU + H/2) + ln0(s) - K^o (]n» + InT). (43) 
2 s + s — 2«7 



Next we deform the integration contour for In 23 in ( f_3"D from a narrow loop around the 
real axis to a wide loop consisting of the two horizontal lines Ims = ±a, with < a < 7. 
The corresponding convolution is denoted by "□" 



0^03 = ^1 Din 03 -In 23, (44) 



and the additional contribution is due to the residue of K\. Taking into account of (|27| , |42| ) 
we find 

lnc = In r- In 03, A In T = A In €, (45) 

where Af(x) = f(x + iO) — f(x — iO) denotes the discontinuity along the real axis. 
Therefore, (f43|) turns into 



J\J n I /j 

In c(s) = — In — + pU + H/2) + In cf>(s) - K X U In 03 - K x o In (_. (46) 

2 s + s — 2z 7 

Lastly, we perform the limit iV — > 00 in the above equations yielding 
lnb = -/?_" + _" 2 n In _ + i_7 o (lnc- In C), 

lnc = -pU/2 + p(fi + H/2) +ln0-^inlnO3-^i'olnC, (47) 
In c = — /5C//2 - /?(/_ + F/2) - In + K\U In 03 + ^ o In (_, 
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where the equation for In c has been derived in analogy to the one for In c, and the function 
has the simplified expression 

In 00) = -20ixy/l - 1/x 2 . (48) 

We want to point out that the function b will be evaluated on the lines Ims = ±a. 
The functions c and c need only be evaluated on the real axis infinitesimally above and 
below the interval [—1, 1]. Also the convolutions involving the "c functions" in fl47D can 
be restricted to a contour surrounding [—1,1] as these functions are analytic outside. 

The detailed derivation of integral expressions for the largest eigenvalue A of the QTM 
is deferred to Appendix A. Here we restrict ourselves to a compilation of the most relevant 
results 




= -2nip-+ [ {In z(s)}' In 1 + < 1 + C ds+ [ [\nz(s - 2i 7 )]' In £(s)ds. 
4 Jc c Jc 



The last two expressions are of particular importance to our further numerical and ana- 
lytical treatment. 

Finally, we want to comment on the structure of the equations determining the thermo- 
dynamical properties of the Hubbard model. In contrast to long-range interaction systems 
[ p| |4"7| ] we have to solve a set of subsidiary equations ([IT]) for the "distribution functions" 
b, c, and c before evaluating the free energy fl49|) . Obviously, the dynamics of the elemen- 
tary excitations of the nearest-neighbor systems is more involved than those of EBL [|7J 
which may be viewed as "free particles with exclusion statistics" . 



6 Numerical Results 



For the numerical treatment of equations (^7J,^) we rewrite them in terms of usual con- 
volutions of functions of a real variable 



K*f 



K{x - y) f(y) dy. 



(50) 



For the functions (27) evaluated on the contours involved in (^TJJ^) we use the notation 
b ± , c 1 * 1 and ^ 



b (x) = b(x ± ia), c (x) = c(x ± i 0), c ± (x) = c(x ± i 0), 



(51) 
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where the shift a in b is arbitrary but fixed with < a < 7. (For many numerical calcu- 
lation we take a = 2^/3 and especially a = 7.) Furthermore, we introduce the following 
relations: 



23± 
Alog£ 



± 



1 + b 
l + c*, 
log(C+/£-), 



23+ 
Alog£ 



l + l/b ± , 

log^/e 7 ). 



(52) 



Thus (J47 ) is written in the form 

log b± = -/3# - K 2i±Q _ Q * log 23+ + K 2i±a+a * log 23' - K 1>±a * A log(c/£), 
log c ± = tff + #i,- Q * log 23+ - Z 1)Q * log 2F + Z lj0 * A log £ ± ± A log £, 
log^± = *f - K h _ a * log 23+ + K 1>a * log 23" - K lfi * A log £ ± ±A log £, 

where 



(53) 



-pU/2 + P(fji + H/2) + log 



±0, 



(54) 
(55) 



< = -PU/2 - /% + - log0 ±o , 
and we have used the notation f a for a function / with shift of the argument by ia 

f a (x) = f{x + ia). 

In particular (j)±Q denotes the function evaluated on the real axis from above/below. 
Notice that the convolution over the terms A log £ and A log £ are determined by Cauchy's 
principal value. Remember that these functions vanish outside the interval [—1, 1]. 
Similarly, from (H9j) we obtain two different relations for the eigenvalue 



log A 



K log[(l 



c+)(l 



c~)] dx 



+ 



-1 
00 



a— 27 



K a ) log 23+ - (/C_ a _ 27 - /C_ a ) log 23-] dx + /3(ji + H/2 + 17/4). 



00 
1 



K log[(l + c + + c+)(l + c~ + c-)/(c+ c-)] dx 



J /C_ 27 log[(l + c+)/(l + c-)]dx -PU/i, 



with 



K{x) 



2tiVT 



X' 



2ixix\/\ — 1/x 2 



(56) 



(57) 



and K, a is the related function with shifted argument. The branch of K, is fixed by the 
requirement K{x) ~ l/(2mx) for large x. By means of the relation 



A log c = -A log + A log £, 
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the first equation of fl5"3|) turns into 

log b ± = V± - K 2 , ±a _ a * log 03+ + K 2 , ±a+a * log QT - if 1)±Q * A log(£/£) , (58) 
where 

*± = -017 - (3H + log ±Q - log ±a _ 27 . (59) 

For the sake of completeness rather than for further applications we mention the results 
for finite Trotter number N. All equations above hold true after the replacement of the 
"driving functions" ip by 



= -m + log ±Q - log ±Q _ 27 - f log 

*± = +(3^ + H/2) + log ±o + f log 



2 to x+so— 217' 



< = -P(ji + H/2) - log ±o + f log 



SQ 



x— so+2i-y ' 



where so is defined in ([H|). These relations for finite Trotter number N have been used 
for a comparison of the results of the integral equations with a direct treatment based on 
the BAE of Sections 3 and 4. Thus it was possible to ensure the accuracy (10~ 6 ) of our 
numerics based on iterations and fast Fourier transform. 

Next we present our numerical results for various physical quantities and discuss them 
in terms of the elementary spin and charge excitations, i.e. "spinons" and "holons" (plus 
gapped excitations based on "doubly occupied sites"). Note that at half-filling the sys- 
tem possesses a charge gap such that the holons do not contribute at low temperatures. 
Furthermore, the hopping integral of the kinetic energy has been set to t — 1. 

In Fig||the temperature dependence of the specific heat is shown for densities n — 1, 
0.8, and 0.5. For half-filling (n = 1.0) the specific heat shows one pronounced temperature 
maximum for lower values of the interaction U. For stronger U this maximum splits into 
a lower and a higher temperature maximum which are due to spin and (gapped) charge 
excitations, respectively. (These findings agree largely with those of f4~8| .) The picture 
remains qualitatively true for small dopings (n = 0.8), however now the lower temperature 
peak receives contributions by gapless charge excitations, hence some weight is shifted 
from higher to lower temperatures. The situation changes quite drastically for fillings 
n m 0.5. Here a pronounced maximum in the specific heat is located at a temperature 
of about T w 0.6 which seems rather insensitive to the interaction. This is explained by 
the irrelevance of the onsite interaction at sufficiently large temperatures, because of the 
low particle density. In addition, we find a maximum at very low temperatures which 
depends very sensitively on U as well as on the particle density n. In order to clarify the 
origin of this additional structure the variation of the specific heat with n is shown in 
Fig.^ for U = 8. Decreasing the particle density n from half-filling (n = 1) to lower values 
(n ~ 0.8) the "spin" maximum at lower temperature increases. This picture is changed 
drastically below n < 0.8. Here the "spin" maximum and its location are suppressed for 
lower n and a shoulder at a higher temperature develops into a clear maximum. This new 
structure in the specific heat is located at about T ~ 0.6 and quite independent of U 
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as already mentioned. We interpret this maximum to be of "charge" type. The complex 
behavior at intermediate densities 0.5 < n < 0.7 is due to a crossover of the "spin" and 
"charge" maxima, see also Fig||. For densities n ~ 1 the "spin" maximum is located at 
finite temperature with finite height whereas the "charge" maximum is located at very 
small temperature with small height. For densities close to n ~ the situation is reversed. 

In Fig.[| and Fig.|6| the magnetic susceptibility x is presented. Again we begin our 
discussion with the half-filled case which is known to correspond to the Heisenberg spin 
chain with interaction strength of order 0(t 2 /U). In fact, we observe a Heisenberg-like 
temperature dependence of the susceptibility with Xmax and T max scaling with U and 
1/U in the range of U = 4,..., 8. Upon doping this behavior remains qualitatively and 
quantitatively unchanged even for n = 0.5. Quite generally, the location T max is shifted to 
lower temperatures, see Fig.|6]. The maximal value Xmax decreases for decreasing particle 
density from n = 1 to n « 0.8. Below the value n < 0.8 the maximum Xmax increases 
for further lowering of the particle density. This behavior is qualitatively explained by 
partially filled bands of charge carriers with spin. 

In contrast to \ the charge susceptibility k (—dn/dfi, i.e. compressibility) shows a 
more interesting dependence on the particle density n, see Fig.|7| and Fig|| At half-filling 
k shows the expected exponentially activated form with vanishing zero temperature value 
due to the charge gap. For any doping the low-temperature behavior is changed completely 
with finite value at zero temperature consistent with a partial filling of the lower Hubbard 
band. For density n = 0.5 we observe two different structures at low temperature similar to 
the case of the specific heat. The lower temperature "spin" peak resembles the structure 
in the susceptibility \i whereas the "charge" maximum at slightly higher temperature 
is caused by the single particle motion of the bare electrons. The compressibility has a 
singular dependence on doping. The smaller the doping the closer the curves are at high 
temperatures and the more divergent at lower temperature, see Fig|| This, of course, is 
exactly the behavior of a system exhibiting a Mott-Hubbard transition at half-filling. 

Our findings are qualitatively in accordance with the results^ of |jTT| , |I~8"|| for the dopings 
treated therein. In particular for specific heats, magnetic and charge susceptibilities the 
results compare well for densities 0.7 < n < 1 and temperatures T > 0.1, giving an 
independent support to the truncation approximation adopted there. 

The present approach has the advantage of explicit evaluations over much wider tem- 
perature and density regions. The T-linearity of the specific heat at very low temperatures, 
as expected from CFT, is clearly observed in contrast to |18|]. Moreover, we have novel 
observations of additional structures at lower temperatures and densities especially in the 
compressibility as mentioned above. These structures can be also interpreted in terms of 
CFT and zero temperature excitations. Therefore, we conclude that the present approach 
is the first one making possible the explicit evaluation of the crossover from the very low 
temperature (CFT) to the very high temperature region in an exact way. 

In Fig.[| we show a separation of the specific heat into spin and charge components. 
This is done in principle on the basis of eigenvalue equations like (OT). As motivated by 



2 In |l7|, [T^] notice the factor 4 for the interaction U. 
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the study of the strong-coupling limit in section 7.1, contributions by b and c functions 
are interpreted as spin and charge contributions, respectively. However, the procedure is 
not unique as we have various alternative formulations resulting in different separations. 
In particular we like to note the expression 

lnA = -/3(e -C/'/4-/i) + J [c A In £/£ - JC ln(l + c+ + c+)(l + c~ + c - )] dx 

/oo roc 
c 2 {x)\n^-(x)dx+ i ci(x)ln^ + (x)rfx, (60) 
-oo J — oo 

with 

Here eo is the groundstate energy at half-filling as given in |] and the additional contribu- 
tions by b and c functions represent correction terms due to spin and charge excitations. 
In Fig.^] we show the results for the specific heat 

where we have applied the separation based on ( |60"D to the temperature derivatives of S and 
n. Note the functional form of the spin part which is rather independent of the doping. 
However, upon small doping the charge contribution develops a low-temperature peak 
which disappears again for larger dopings. We would like to comment on the issue of the 
"mathematical separation" of spin and charge as described above (and similarly applied 
in [|17[ [Tj|]) that it may give rise to artificial results. For instance, at higher temperatures 
the "partial specific heats" show negative values whereas the total specific heat is always 
positive. In section [T3] the spin-charge separation is treated at low temperatures and 
arbitrary particle density via an involved interplay of the various degrees of freedom 
rather than by a superficial interpretation of formulas. 



7 Analytical solutions of the integral equations 

In the previous sections we have derived non-linear integral equations for the largest eigen- 
value being directly related to the free energy of the Hubbard model at finite temperatures 
T = 1/(3. For arbitrary temperatures and densities the integral equations can be solved 
only numerically. However, in some limiting cases analytical results can be derived and 
relations obtained which permit a comparison to known analytical results. This implies 
the consistency of our approach. 
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7.1 Strong-coupling limit 

In the strong-coupling limit U — > oo at half-filling (// = 0) the Hubbard model is expected 
to reduce to the Heisenberg chain. Indeed, in the strong-coupling limit we find directly 
the thermodynamics of the Heisenberg model. This can be seen as follows: 

Considering the limit 7 — >• 00 (with 7 = {7/4) we rescale the argument of the auxiliary 
functions by x h> 27a; as well as the ratio fl/(2j) 1— > (3 and (2 r y)H 1— ► H . It turns out that 
all contributions of c 1 * 1 and can be dropped in ( |58D because of their vanishing range of 
integration. Moreover c ± and c ± tend to zero. The remaining equations read 

log b ± = ^- K 2)±a . a * log <B + + K 2 ,± a+a * log 9T with ± = -PH + 2<K~pK x , ±a . 

(63) 

According to equation fl56|) and after dropping the groundstate energy shift /3{7/4 the 
eigenvalue is 

logA = K 1} _ a * log<B + U =0 - K 1>a * log<B"U =0 + /3i7/2. (64) 

Now we define 

b:=l/b + , B:=l + b and b := b~, S := 03", (65) 
which are inserted into the integral equations. By means of the identity 

log 93 + = log B — log 6, 

we have 

- log b - K 2fi * log b = +\E^" - K 2fi * log B + if 2 ,2a * log -B, 

- log b - K 2fi * log b = -\I> b + K 2 -2a * log S - K 2 ,o * log 5. 

Using the Fourier transform this provides 

log b = +/3H/2 - 2vr/3 $ +Q + R * log B - J7 +2a * log 7J, (66) 
log 6 = -{3H /2 - 2ir{3 <$>^ a + Ro* log B-R_ 2a * log B, (67) 

with the eigenvalue 

/oo 
( $ a log £-$_ Q logS) dx, (68) 
-00 



and 



= ; r^r, Rn 



dk e 



ikx—ka 



2 sinli7r(a; + ia) ' a y_ oc 27T 1 H- e' fc l 

These relations correspond to the non-linear integral equations of the isotropic antiferro- 
magnetic spin- 1/2 Heisenberg chain [p5, |26|. In addition, this case is also related to the 



thermodynamics of the t — J model at half-filling [30 



19 



7.2 Free-Fermion limit 

Let us consider the opposite limit U — > (at \i — and H = 0), that is, the case of two 
independent Free-Fermion systems. The non-linear integral equations (^) simplify to an 
algebraic set of equations 

log b ± = -(3H - log <B + + log <B~ - (| ± ~) A log(c/€) , 
log c± = + if/2) + log 0± - log ®~ + (| ± |) A log C, 
logF = -/3(// + H/2) - log ± - log 23+ + (-i ± |) A log €, 



with 



, . i' 1 k- , (l + c+ + c+)(l + c- + c-)(l + , 
log A = - / /Co log ^ =^=- — ^ dx. 

-1 C+C"(l + C+) 



The solution reads as follows 



e (3H /(f) _|_ e 2/3 M + J + e /3( M -H/2) ' 



C+ 



1 _|_ ePQi-H/2) /0 i + b+ ' 1 + b~ 

1 l+e^~ H / 2 )/0 1 



e /3(^+Jf/2)0^ + ' e P(K.+H/2) 1 _|_ e /3(M-H/2)0 1 _j_ ' 

Lastly, we substitute x = sin in the integration for the eigenvalue which leads to 

1 /' 7r 

logA = + — / ln[l+exp(/3(// + .ff/2 + 2cosfc))] dfc 

1 P 

+ — / In [1 + exp (p(ji - H/2 + 2 cos k))} dk. 

2vr ./_„ 

This is the desired result. 



7.3 Low-temperature asymptotics 

The low-temperature regime is the most interesting limit as the system shows Tomonaga- 
Luttinger liquid behavior. We will derive analytic expressions for the thermodynamics 
within our first principles calculations and confirm the field theoretical predictions. In 
particular we will show how the non-linear integral equations correspond to the known 
dressed energy formalism of the Hubbard model. This represents a further and in fact the 
most interesting consistency check. 
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For T = < we can simplify the non-linear integral equations as follows. We 
adopt fields H > 0, \i < 0, such that b — — ^ 0, c^ 3 — > and l/c~ — > at (3 3> 1 which can 
be verified numerically. Moreover, one finds 

|b + |, |l/c + | 3> 1 for \x\ < A , o"o and |b + |, |l/c + | <C 1 for |x| > A , cr , 

for certain crossover values Ao, cro- The slopes for the crossover are steep, so that the 
following approximations to the integral equations (|53| ) are valid: 

•rXo r+h) 



log b+ = <f) b - f ° K 2 {\ - A') log b+(A') dA' + / ° K ha {\ - sin k!) cos k' log c v (A;') dk' 

/+Ao 
#i _ Q (sin k - A') log b + (A') dA'. 
-An 



(69) 



Here we use c v (k) = l/c + (sin/c) and <To = sinfc . The driving terms read 

7T 2 (ir 2 (A-Ao)+K 2 (A + A )) 



fe (A) = -(3e° s (\) - 
+ 



6(logl/b+)'(A ) 
7r 2 cos k (K la (\ — sin fc ) + K la (\ + sin fc )) 
6(logl/c v )'(A; ) • 



a/W- a e 0( \\ , 7r 2 (^i,-q(smfc - A ) + i^ 1; _ a (sinfc + A )) 

C (A) - -/fe e (A) + 6 (logl/b+)'(A ) ' (70) 

where e® = H, e° c = —\i — U/2 — H/2 — 2cos/c. Retaining only the leading terms in the 
integral equations and choosing the imaginary part of the integration contour as a = 7, 
we find the following connections between auxiliary functions and the dressed energy 
functions: 

logb+ = -pe, + 0(1/(3) and log c v = -(3e c + 0(1/(3). (71) 



For a comparison with fEq , [13[ note the different normalization of the chemical potential. 
The free energy also admits the same approximation scheme yielding up to 0(T 2 ) terms 
in the low-temperature expansion. To present this, we introduce "root density functions" 
p by the definition 



K 2 (\ - A') Ps (\>) dA' + / K ha (X - sin k') p c (k') dk' 

-Ao J—ko 



+a ( 72 ) 



^ /-Mo 

pc(k) = hcosA; / K\ _ a (sin k — A') p s (A') dA'. 

27r J-* 



Note that the kernel matrices for the integral equations (|)9|, [72|) are mutually transpose. 
The following equality is an immediate consequence: 
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1 

2^ 



+h> r+ko /"+Ao 

\ogc v (k)dk = Pc (k)<f )c (k)dk+ Ps (X)<j )b (X)dX. (73) 

fco J—ko J— Ao 



With these relations the eigenvalue (log A) reads 

»+fco /"+Ao 

p c (fc) c (fc) dA; + / ,,. 

-fc ./-Ao 6^(A;q 



logA = -/3 7 +/ Pc (k)<p c (k)dk+ I p s (\)fo(\)d\+ —?- , (74) 



where we have replaced (log l/c v )'(/c ) in the denominator by /3^(fc ). Substituting Q70D 
into ([74]) we arrive at the final expression, 

/ = -o-^(- + -V (75) 



The definitions of sound velocity and the ground state energy coincide with standard 
results, v S)C = £' SjC /2irp s , c \x ,k , and 

/+ko /"+Ao 
Pc (k)e° c (k)dk+ / p s (X)e° s (X)dX. 
■fco ^-A 

Here the trivial shift in the energy U/4 is omitted. We thereby conclude that our formalism 
completely recovers the correct contribution from spinon and holon excitations in the low 
temperature behavior. This is a manifestation of spin-charge separation due to which 
each elementary excitation contributes independently to (|75| ) where the velocities v c and 
v s typically take different values. 

7.4 High-temperature limit 

Finally, we consider the high-temperature limit T — > oo with H, U as well as ftp, fixed. 
The auxiliary functions in fl58|) become constant 



logb ± = 0, log c ± = f3 p — log 2, log c ± = — {3 p — log 2, 
log A = log(l + c+ c)/c. 



Thus, the free energy reads 



/= -2Tlog(l + e M/T ) with p/T = log-^—, (76) 

2 — n 

where n is the particle density. We obtain the entropy 

2 n 

S = 2\og- nlog- , (77) 

2 — n 2 — n 

as expected by counting the degrees of freedom per lattice site. Especially, at half-filling 
n — 1 this equals to S = log(4). 
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8 Summary and Discussion 



In this paper, the novel formulation of thermodynamics for ID quantum systems has been 
successfully applied to the Hubbard model. Several quantities of physical interests have 
been evaluated with high numerical precision and various limiting cases have been studied 
analytically. 

As already noted above, we may consider as one of the most practical advantages of 
the present formulation the fact that one only has to deal with a finite number of un- 
known functions and nonlinear integral equations among them. This does not only imply 
convenience, rather it opens a more fundamental understanding related to the particle 
picture of ID quantum systems. For the Heisenberg model, the complex conjugate auxil- 
iary functions play a role which seems to correspond to the elementary spinon excitations 



50| 51], while for the integrable t — J model one further function related to the holon 
is needed. In this paper, we have shown that three independent functions b, c, c describe 
the complete thermodynamics, physically corresponding to spinons, and holons in upper 
and lower Hubbard bands. In the T — > limit, these functions are shown to reduce to 
energy density functions ("dressed energy functions") for such elementary excitations. 

This interpretation which is natural at low temperatures poses however a problem at 
finite temperatures. The auxiliary functions, related to energies of excitations at T = 0, 
are no longer real for arbitrary temperature. Thus they lose the direct connection to 
physical excitations in the sense of energy levels. On the other hand, imaginary parts of 
energies indicate a finite life-time of excitations, or in this case a decay of the elementary 
particles of the system. We leave the investigation of these questions as an interesting 
future problem. 

Obviously, our formulation can be extended to the evaluation of the asymptotics of 
correlation functions, such as spin-spin correlation lengths etc. These investigations on 
highly correlated electron systems, including the Hubbard, the integrable t — J, super- 
symmetric U models will be reported together in the near future. 
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A Derivation of integral expressions for the eigen- 
value 



Here we turn to the derivation of expressions for the largest eigenvalue of the QTM 
in terms of the above auxiliary functions. We calculate Y2j ^ n z j by a Cauchy integral 
of the function f(s) = lnz(s) [In (1 + h/h(s))]' where the zeros of (1 + h/fa(s)) in the 
neighborhood of the real axis are precisely the Sj. Furthermore we use a contour £ 
surrounding the Sj in anticlockwise manner. The Sj are not located on the branch cut 
of lnz(s) from —1 to 1, hence £o consists of two disconnected parts. (For not too low 
temperatures and vanishing external fields these parts are loops around ] — oo, —1] and 
[l,oo[, respectively. In the general case they are appropriately deformed.) However the 
Zj corresponding to a particular Sj might have to be calculated by the use of the second 
branch of In z(s). We therefore obtain 



2mJ2^z( Sj )= f(s)\ lhianch ds+ f(s)\ 2hr&nch ds, 

JCo JCo 



(7c 



where the first and second term on the right hand side will be abbreviated by £i and S 2 , 
respectively. We next manipulate Si 



Si 



In z(s) 



Co=-(Ci+C2+C 3 



ds, 



In z(s) 



Ci+C 3 



In ( 1 + |( a ; 



ds + J In z(s) [hi T(s)] ds, 



(79) 



where in the first line the integration contour £ can be replaced due to Cauchy's theorem 
by three contours (taken in anticlockwise manner): L\ from zoo to —1, surrounding [—1,1], 
from —1 back to 200; £2 surrounding the axis Im(s) = —7 (where the simple poles of 
[1 + U/ls(s)] are located, which are identical to the simple zeros of T(s)); £3 around so 
(which is a pole of order N/2 of [1 + U/h(s)]), see (|3T|) . Next we replace £2 by contours 
£1, £— 2«7 (where £ surrounds the real axis), and £3—227 



£2 



\nz(s) [in T(s)] ds 



\nz(s) [\nT(s)]'ds 

/ In z(s - 2i 7 ) [In T{s - 2iy)] ' ds. 
Jc 3 +c 



10) 



In the last term of (|S0|) the function lnT(s — 2ij) can be replaced by — ln*B(s) as the 
difference of these functions amounts to an analytic contribution vanishing in the contour 
integration. Of course, the £3 integration in (|T5D and ( j5D| ) can be done explicitly yielding 



Si = niN \n[z(so)z(so — 2^7)] 
In z(s 



in (t(s) (1 + |( S ; 



ds + J lnz(s-2i~f) [In <B(s)]' ds. 
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Next we perform integration by parts on the right hand side where the first integral also 
contributes a "surface term" 



Ei = -2ttz In [(1 + e -^ +H ^) (1 + e ^~ H / 2 ))] + iriN \n[z(s )z(s - 2ij)] 



+ [ [lnz(s)]'\n 



tw(i + ^ : 



ds - J [\nz(s - 2ij)]'\n<B(s)ds, 



f82 x 



where now the integration along £1 can be restricted to a loop £ 4 along the cut from -1 to 
1. In the limit N — > oo we can replace 7ix/Vln[z(so)£(so — 2*7)] by 27riA^ln + 47ri/?7 

such that in combination with (|24l) 



2vrilnA = 2vri/?(/i + H/2 + U/A) + S 2 + / [lnz(s)]'ln 

J£ 4 



( 1 + 



-y pnz(s-2i7)]'ln?8(s)ds. 

A very similar line of reasoning yields 

2m]nA = 2m0(n- H/2 + U/A) + S' 2 + f [In In T(s) ( 1 + 

-y [In z(s + 217)]' In »(s)ds, 

where S' 2 is defined similarly to £2 after interchanging Z3 and I4. Combining 
find the symmetrised version 

Am In A =4vrz/5(/i + U/A) + S 2 + E 2 



33) 



34) 



we 



+ / [ln^(s)]'ln 
J U 



T(,)T(8) ( 1 + 



1 + ^ 



*5) 



~y [In 2(5-2^7)]' In <B(s)ds-y [lnz(s + 217)]' In »(s)ds. 



The present formula is still inconvenient as the first terms on the right hand side contain 
"non-standard" functions. However, we can substitute the terms 



1 + {*<»: 



14 



1 + ^: 



1+^: 



(86) 



2. branch 



1. branch 

without change of the integral as [In z(s)]'|i. branch = — [In z(s)]'|2.branch- For the same reason 
we find 



[lnz(s)]'ln 



1 + 



2. branch 



57) 



2. branch 



25 



where functions have to be evaluated on the 1. branch unless indicated differently. Insert- 
ing ( |36"l , j37|) into (jS5p, combining the contours £ an d £4 info C, and simply extending £ 4 
to C for the integrals involving T and T (due to analyticity) we arrive at 



A-ni In A =4?ri/3(// + U /A) 
+ / [lnz(s)]'ln 



T{s)T{s) 



{4 



2. branch- 

- / pnz(a- 2*7)]' In »(s)ds- / [In z(s + 2ij)]' In lB(s)ds. 



Next we find the identity 



(U/h 



I 2. branch 



(bT/T) | x 



branch 



Hence the integrand of the first integral in is 



TT 



bT\ ( T 

1+ t) [ 1 + w 



(t + bry 



[1 + c + cr®®, 



(89) 



(90) 



with all functions on the first branch. We are now in the position to formulate the first 
main expression for the eigenvalue 



Am In A = Aixi(3{{i + U/A) + 2 J [In z(s)}' In (1 + c + c) ds 



In 



z(s — 2ry) 



z(s 



ln03(s)ds 



In 



z(s + 2ry) 



z(s 



(91) 



\n<8(s)ds. 



The last equation can be reduced further by substituting the 03 integral by a 03 integral. 
For this purpose we employ 



In 



z(s + 2ij) 

z(s + 227) 



z{s) 
z(s + 2ij) 



ln03(s)ds 

ln(h + l 2 + h + h + h + k + h + h)ds 
hi(r i +T2 + h+T A )ds. 



(92) 



In the first term on the right hand side the contribution due to z(s + 227) vanishes upon 
taking the contour integral and can be replaced by (the equally vanishing) z(s — 227). In 
the second term the contour is replaced by —(C—2i^). Using (Zi + h + h + h){s — 2ry) = 
exp(-pH)(h + h + h + U)(s)/[<j) + (s)<f)-(s)] we find 



L 



L=-{L-2i-i) 



hi 



z(s + 227) 



Z(S 



\n(h + k + k + h)ds 

z(s — 2ry) 



2nipH + 



In 



z(s 



(93) 



ln(/i + l 2 + l 3 + h)ds, 
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where terms involving + have been dropped as they vanish in the limit N 
Inserting this into (|92|) we obtain 



oo. 



In 



z(s + 2ij) 

4s) 



ln<B(s)ds = -2nif3H + 



In 



z(s — 2ij) 



Z(S 



ln»(s)ds, 



(94) 



and with (|91~D we find 



2ttz In A = 2rrip(fi + H/2 + U/A) + j [In z(s)]' In (1 + c + c) ds 

z(s — 2ij) 



In 



Z(S 



(95) 



ln<B(s)ds. 



Finally, we want to show how to express the eigenvalue entirely in terms of c functions, 
i.e. without contributions by b. To this end we note the identity 



[In z(s)}' \nr{s)ds = - 2m(3{fi + H/2 + U/2) 



£=-(£-217) 



— J [In z(s — 2ry)]' In r(s — 2ij)ds, 



(96) 



where we have dropped terms that do not contribute in the limit N — > oo. Next we replace 
the t(s) and r(s — 2ry) functions by 



T = m, [t(s - 2?7)] _1 = r(s) 



/x + / 2 



^3 + ^4 + Zi + / 2 + + h 



(97) 



Since the function (/ x + l2)/(h + h + h + h + h + U) is analytic in the neighborhood of 
the real axis its contribution to (|96h vanishes. Therefore (PB?) results into 



J [In z(s))' (In c(s) + In <B(s))ds = - 2vr^(/i + + C//2) 



+ J [\nz(s -2i 7 )]'(rn<8(s) +ln£(s))ds. 



(98) 



Inserting this into (|95|) we are left with 



27TzlnA = -2mP T + [ [In z(s)}' In 1 + ^ + C ds+ / [In z(s - 2i 7 )]' In £(s)ds. (99) 
4 ./ /• c 
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Figure 3: Specific heat versus T for particle densities n — 1, n — 0.8 and n = 0.5. 




Figure 4: Specific heat versus T for fixed U = 8. 
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n^l.OO 



n=0.80 



n=0.50 




Figure 5: Magnetic susceptibility versus T for n — 1, n — 0.8 and n = 0.5. 




Figure 6: Magnetic susceptibility versus T for fixed U — 8. 
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Figure 7: Compressibility versus T for particle densities n — 1, n — 0.8 and n = 0.5. 




Figure 8: Compressibility versus T for fixed U — 8. 
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Figure 9: Separation of specific heat (solid) in spin (dashed) and charge components 
(dotted). 
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